\defmodule {GammaProcess}

This class represents a \emph{gamma} process \cite[page 82]{fMAD98a}
$\{ S(t) = G(t; \mu, \nu) : t \geq 0 \}$ with mean parameter $\mu$ and
variance parameter $\nu$.
It is a continuous-time process with stationary,
independent gamma increments such that for any $\Delta t > 0$,
\begin{equation}
 S(t + \Delta t) = S(t) + X,\label{GammaEqn}
\end{equation}
where $X$ is a random variate from the gamma distribution
\texttt{Gamma}$(\mu^{2}\Delta t / \nu , \mu / \nu)$.

In this class, the gamma process is sampled sequentially using
equation (\ref{GammaEqn}).

\bigskip\hrule\bigskip

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{code}
\begin{hide}
/*
 * Class:        GammaProcess
 * Description:
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       Pierre Tremblay, Jean-Sebastien Parent & Maxime Dion
 * @since        July 2003

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
\end{hide}
package umontreal.iro.lecuyer.stochprocess;\begin{hide}
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.randvar.*;
import umontreal.iro.lecuyer.probdist.GammaDist;
import umontreal.iro.lecuyer.util.Num;
\end{hide}

public class GammaProcess extends StochasticProcess \begin{hide} {
   // Make sure that process is strictly increasing: when two successive steps
   // are equal, reset the new one to   old*(1 + EPS)
    protected static final double EPS = 1.0e-15;
    protected boolean      usesAnti = false;

    protected RandomStream stream; // used in the generatePath(), so must be kept
    protected GammaGen     Ggen;
    protected double       mu,
                           nu;
    protected double       muOverNu,
                           mu2OverNu;

    // For precomputations for standard GP
    protected double[]     mu2dtOverNu;

   protected void setLarger (double[] path, int left, int mid, int right) {
      /* make sure that path[left] < path[mid] < path[right] by adding a
         small epsilon to the last two */
      double myEps;
      if (path[left] >= 0)
         myEps = EPS;
      else
         myEps = -EPS;
      path[mid] = path[left] * (1.0 + myEps);
      if (path[mid] >= path[right])
         path[right] = path[mid] * (1.0 + myEps);
   }


   protected double setLarger (double[] path, int left, int right) {
      /* make sure that path[left] < s < path[right] by adding a
        small epsilon to the last two; return s. */
      double myEps;
      if (path[left] >= 0)
         myEps = EPS;
      else
         myEps = -EPS;
      double s = path[left] * (1.0 + myEps);
      if (s >= path[right])
         path[right] = s * (1.0 + myEps);
      return s;
   }


   protected double setLarger (double v) {
      /* return a number that is slightly larger than v */
      if (v >= 0)
         return v*(1.0 + EPS) + Num.DBL_EPSILON;
      else
         return v*(1.0 - EPS);
    }\end{hide}
\end{code}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Constructors}
\begin{code}

   public GammaProcess (double s0, double mu, double nu,
                        RandomStream stream)\begin{hide} {
        this (s0, mu, nu,  new GammaGen (stream, 1.0));
// Note that the parameters (the 1.0) supplied to
// the GammaGen constructors are meaningless.
// The 'real' parameters are supplied at time of generation
    }\end{hide}
\end{code}
\begin{tabb} Constructs a new \texttt{GammaProcess} with parameters $\mu =
\texttt{mu}$, $\nu = \texttt{nu}$ and initial value $S(t_{0}) = \texttt{s0}$.
The gamma variates $X$ in (\ref{GammaEqn}) are generated by inversion
using \texttt{stream}.
\end{tabb}
\begin{code}

   public GammaProcess (double s0, double mu, double nu, GammaGen Ggen) \begin{hide} {
        this.mu    = mu;
        this.nu    = nu;
        this.x0    = s0;
        this.Ggen     = Ggen;
        stream = Ggen.getStream();
    }\end{hide}
\end{code}
\begin{tabb} Constructs a new \texttt{GammaProcess} with parameters $\mu =
\texttt{mu}$, $\nu = \texttt{nu}$ and initial value $S(t_{0}) = \texttt{s0}$.
The gamma variates $X$ in (\ref{GammaEqn}) are supplied by the gamma random
variate generator \texttt{Ggen}. Note that the parameters of the
\externalclass{umontreal.iro.lecuyer.randvar}{GammaGen} object
\texttt{Ggen} are not important since the implementation forces the generator
to use the correct parameters (as defined above).
\end{tabb}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Methods}
\begin{code}\begin{hide}
   public double nextObservation()  {
        double s = path[observationIndex];
        double v = s;
        s += Ggen.nextDouble (stream, mu2dtOverNu[observationIndex], muOverNu);
        if (s <= v)
             s = setLarger (v);
        observationIndex++;
        observationCounter = observationIndex;
        path[observationIndex] = s;
        return s;
    } \end{hide}

   public double nextObservation (double nextT) \begin{hide} {
        // This method is useful for generating variance gamma processes
        double s = path[observationIndex];
        double v = s;
        double previousT = t[observationIndex];
        observationIndex++;
        observationCounter = observationIndex;
        t[observationIndex] = nextT;
        double dt = nextT - previousT;
        s += Ggen.nextDouble (stream, mu2OverNu*dt, muOverNu);
        if (s <= v)
             s = setLarger (v);
        path[observationIndex] = s;
        return s;
    }\end{hide}
\end{code}
\begin{tabb} Generates and returns the next observation at time $t_{j+1} =
\texttt{nextTime}$,
using the previous observation time $t_{j}$ defined earlier
(either by this method or by \texttt{setObservationTimes}),
as well as the value of the previous observation $X(t_j)$.
\emph{Warning}: This method will reset the observations time $t_{j+1}$
for this process to \texttt{nextT}. The user must make sure that
the $t_{j+1}$ supplied is $\geq t_{j}$.
\end{tabb}
\begin{code}

   public double[] generatePath() \begin{hide}  {
        double s = x0;
        double v;
        for (int i = 0; i < d; i++) {
            v = s;
            s += Ggen.nextDouble(stream, mu2dtOverNu[i], muOverNu);
            if (v <= s)
               s = setLarger (v);
            path[i+1] = s;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    } \end{hide}
\end{code}
\begin{tabb}
Generates, returns and saves the path
$\{X(t_{0}), X(t_{1}), \ldots, X(t_{d})\}$.
The gamma variates $X$ in (\ref{GammaEqn}) are generated using the
\externalclass{umontreal.iro.lecuyer.rng}{RandomStream} \texttt{stream}
or the \externalclass{umontreal.iro.lecuyer.rng}{RandomStream} included
in the \externalclass{umontreal.iro.lecuyer.randvar}{GammaGen} \texttt{Ggen}.
\end{tabb}
\begin{code}

   public double[] generatePath (double[] uniform01) \begin{hide} {
        double s = x0;
        double v;
        for (int i = 0; i < d; i++) {
            v = s;
            s += GammaDist.inverseF(mu2dtOverNu[i], muOverNu, 10, uniform01[i]);
            if (v <= s)
                s = setLarger (v);
            path[i+1] = s;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    } \end{hide}
\end{code}
\begin{tabb} Generates, returns and saves the path $
\{X(t_{0}), X(t_{1}), \ldots, X(t_{d})\}$. This method does not use the
\externalclass{umontreal.iro.lecuyer.rng}{RandomStream} \texttt{stream} nor the
\externalclass{umontreal.iro.lecuyer.randvar}{GammaGen} \texttt{Ggen}. It
uses the vector of uniform random numbers $U(0, 1)$ provided by the user
and generates the path by inversion.  The vector \texttt{uniform01} must be of
dimension $d$.
\end{tabb}
\begin{code}

   public void setParams (double s0, double mu, double nu) \begin{hide} {
        this.x0 = s0;
        this.mu = mu;
        this.nu = nu;
        if (observationTimesSet) init(); // Otherwise no need to.
}\end{hide}
\end{code}
\begin{tabb} Sets the parameters  $S(t_{0}) = \texttt{s0}$,
$\mu = \texttt{mu}$ and $\nu = \texttt{nu}$ of the process.
\emph{Warning}: This method will recompute some quantities stored internally,
which may be slow if called repeatedly.
\end{tabb}
\begin{code}

   public double getMu() \begin{hide} { return mu; }\end{hide}
\end{code}
\begin{tabb}
Returns the value of the parameter $\mu$.
\end{tabb}
\begin{code}

   public double getNu() \begin{hide} { return nu; }\end{hide}
\end{code}
\begin{tabb}
Returns the value of the parameter $\nu$.
\end{tabb}
\begin{code}

   public void setStream (RandomStream stream) \begin{hide} {
        this.stream = stream;
        Ggen.setStream (stream);
    }\end{hide}
\end{code}
\begin{tabb}
Resets the \externalclass{umontreal.iro.lecuyer.rng}{RandomStream}
of the \externalclass{umontreal.iro.lecuyer.randvar}{GammaGen}
to \texttt{stream}.
\end{tabb}
\begin{code}

   public RandomStream getStream() \begin{hide} { return stream; }\end{hide}
\end{code}
\begin{tabb}
Returns the \externalclass{umontreal.iro.lecuyer.rng}{RandomStream} \texttt{stream}.
\end{tabb}
\begin{code} \begin{hide}

    protected void init() {
    // init of StochasticProcess only sets path[0] = x0 if observation times are set.
        super.init();
        muOverNu  = mu / nu;
        mu2OverNu = mu * mu / nu;
        mu2dtOverNu = new double[d];
        if (observationTimesSet) {
            for (int i = 0; i < d; i++)
                mu2dtOverNu[i] = mu2OverNu * (t[i+1] - t[i]);
        }
    }

}\end{hide}
\end{code}
